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Understanding the conditions ensuring the persistence of a population is an issue of primary im¬ 
portance in population biology. The first theoretical approach to the problem dates back to the 
50’s with the KiSS (after Kierstead, Slobodkin and Skellam) model, namely a continuous reaction- 
diffusion equation for a population growing on a patch of finite size L surrounded by a deadly 
environment with infinite mortality - i.e. an oasis in a desert. The main outcome of the model 
is that only patches above a critical size allow for population persistence. Here, we introduce an 
individual-based analogue of the KiSS model to investigate the effects of discreteness and demo¬ 
graphic stochasticity. In particular, we study the average time to extinction both above and below 
the critical patch size of the continuous model and investigate the quasi-stationary distribution of 
the number of individuals for patch sizes above the critical threshold. 

PACS numbers: 


I. INTRODUCTION 


Many biological or chemical processes involve the dy¬ 
namics of discrete “particles” (e.g., molecules or organ¬ 
isms) that diffuse and interact with each other and/or 
with an external environment [IH3. In population dy¬ 
namics, for instance, individuals spread in space, inter¬ 
act, reproduce and die. 

If the number density of particles is very large, the 
macroscopic description in terms of continuous fields is 
typically appropriate. A well established approach to 
model the spatio-temporal evolution of the population 
density field, V{x,t), is in terms of a reaction-diffusion 
equation that, in one spatial dimension, reads 


^ - B — 
dt dx^ 


+ fir ). 


( 1 ) 


In the above equation, D is the diffusion coefficient and 
/{V) rules the chemical kinetics or, in the language of 
population biology, the local rate of population growth. 
A well known growth term is the logistic model 


m = rv (i-|) , 


( 2 ) 
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whose quadratic term models the competition for re- 
sources, K being the local carrying capacity (i.e. the 
maximal population density that can be sustained) while 
r denotes the intrinsic growth rate. 

Conversely, if the density of particles is not very large 
the discrete nature of the population cannot be ne¬ 
glected Q, new effects arise and the continuous descrip¬ 
tion becomes inaccurate. In order to account for the new 
effects several approaches are possible. For instance, for a 
population of N individuals, a possibility is to introduce 
a cutoff at the value 1/N of the normalized density for 
the continuous field equations 0 (see for a review). 
Another possibility is to supplement Eq. m with a noise 
term accounting for microscopic fluctuations originated 
by the finite number of individuals dMl. A further 
modeling step is to consider a lattice where each site is 
occupied by an integer number of individuals [l^ . [l5l| . or a 
contact process (see, e.g., dill and references therein), 
where each site is occupied at most by an individual. 

In the present work, we consider a system of particles 
diffusing in space and interacting when they get within 
a given interaction distance. In particular, revisiting the 
continuous KiSS model, after Kierstead, Slobodkin and 
Skellam dEl, we investigate the effects induced by 
discreteness on the dynamics of a population inhabiting 
a favorable region (“an oasis”) surrounded by a deadly 
environment (“a desert”). Accounting for such effects is 
important to assess the role of demographic and environ¬ 
mental stochasticity on extinction dynamics 
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Several efforts have already been undertaken along this 
direction within the framework of the KiSS model. Some 
model modifications that include effects of intrinsic vari¬ 
ability on the critical habitat size ensuring population 
persistence have been studied in relation to environmen¬ 
tal stochasticity [2l|. As for demographic stochasticity 
due to population discreteness it has been investigated 
either using a cutoff on the continuous field equations 
( 2 ^, in the same spirit of Ref. Q, or within a stochas¬ 
tic framework using field theory techniques [ 2 ^. In our 
work, we adopt a discrete particle model [IJ focusing on 
demographic stochasticity effects. Moreover, our study 
is not restricted to steady state properties but also deals 
with dynamical features of extinction phenomena. 

The article is organized as follows. In Section [IT] we re¬ 
call the continuous KiSS model [ll,[l3 and introduce its 
discrete analogue. In Section IIIII we numerically study 
the long time properties of the discrete model, charac¬ 
terizing the so-called quasi-stationary state [1^, , and 

comparing it with the properties of the continuous model. 
Section |IV] is devoted to the problem of determining the 
critical habitat size for population persistence and to the 
estimation of extinction times in the discrete particle sys¬ 
tem. Discussions and conclusions are presented in Sec. El 
As a side discussion, in the Appendix we briefly consider 
the effects of changes of the microscopic rules of the dis¬ 
crete model. 


II. MODEL 

A. Continuous KiSS model 

Determining the conditions for persistence or, con¬ 
versely, estimating the probability of extinction of species 
and populations is a question of paramount importance 
in population biology. The theoretical treatment of the 
problem was initiated by Kierstead, Skellam and Slobod- 
kin [ 1 ^, [ 1 ^ who considered a population that, ruled by 
Eqs. (IT][5]), grows with rate r and diffuses with diffusion 
coefficient D within a patch of size L surrounded by a 
deadly environment. For a given growth rate, higher dif- 
fusivities imply larger fluxes across the boundaries, so 
that larger patches are necessary to compensate the pop¬ 
ulation loss in order to allow stable persistence. It is 
then natural to determine the minimal (critical) patch 
size, Lc, ensuring population persistence. At least di¬ 
mensionally, the critical patch size can be obtained by 
balancing the growth rate r with the exit rate from the 
patch due to diffusion This way one derives 

that Lc ^ \/D/r that is the correct result up to an order 
one constant (see below). 

In one-dimensional space, for simplicity, it is enough 
to consider Eqs. (fT1l2l) in the interval 0 < x < L, rep¬ 
resenting the favorable patch, while the outside hostile 
environment results in the boundary conditions 


(3) 


Detailed analysis of the KiSS model can be found in 
[TsL [T^ . , for generalizations of the model includ¬ 

ing different kind of population growth terms see pin. 
Here, we only recall the main results on the critical patch 
size. We notice that, close to the extinction, the popu¬ 
lation density will be very small {V <C K) so that we 
can linearize the growth term (2)) by posing f{V) = rV. 
Starting from a generic initial condition, population den¬ 
sity will thus grow (decay) depending on the positive 
(negative) sign of the largest eigenvalue Ai of the (lin¬ 
ear) operator C = Dd1 + r with boundary conditions ([2]). 
Standard computation shows that the eigenvalues are 

n = l,2,... (4) 



with Lc = TT\/D/r. Consequently, the population can 
persist only if the favorable habitat is larger than the crit¬ 
ical patch size Lc- However, while in the linear case also 
the eigenfunctions and thus the stationary population 
density can be easily derived, for the logistic case only 
approximate solutions are known (see, e.g.. Ref. M)- 

It is worth remarking that proper estimations of Lc, 
with suitable generalizations of the above model, are used 
in the design of protected areas for endangered species, 
in the context of conservation biology [33 - l33 | . 

In the following we focus on the logistic growth case. 
For later convenience it is useful to formulate the contin¬ 
uous model (mi2]) in non-dimensional variables. From the 
above analysis it is natural to measure time in units of 
the inverse growth rate 1/r and space in units of \JD jr. 
Finally, normalizing the population density to the carry¬ 
ing capacity, that is introducing 9 = VlK^ we can rewrite 
Eqs. (Hill]) as 


dt 




(5) 


with 0(0) = 0(L) = 0, where L now is the non- 
dimensional domain size. With these units the critical 
patch size is Lc = tt. 


B. Discrete KiSS model 

In recent years, several efforts have been devoted to 
the extinction problem in populations with a finite num¬ 
ber, N, of individuals, using non-spatial models 
[ 33 . The factors influencing extinction, through fluctu¬ 
ations and decline of a population, are in principle ex¬ 
tremely varied and attributable to a wide class of biologi¬ 
cal and environmental causes, among which demographic 
stochasticity [s^], associated to the unavoidable random 
occurrence of birth and death events, is one of the most 
important ones. 

In the context of spatially extended systems, some 
effects of demographic stochasticity on steady patterns 
of populations have been analyzed using different mod¬ 
els [12,[ 2 ^. In particular, in Ref. the discreteness was 


r{0)=r{L) = 0. 







3 


modeled using reaction-diffusion equations with a cutoff 
on the growth term, and it was found that the critical 
patch size for survival is larger than in the model with¬ 
out cutoff. A similar effect was found also in Ref. [^ . 
using a stochastic discrete model, where it is shown that 
the population can go extinct above the critical patch 
size of the continuous model and an “effective” (finite) 
critical patch size can be introduced only when the com¬ 
petition term is weak enough. In our study, similarly to 
, the discrete population can go extinct for habi¬ 
tat sizes guaranteeing the persistence in the continuous 
case. However, we will always refer to the critical patch 
size of the continuous model Lc because, as we will see, 
this value still marks a transition in the system behavior. 

We now introduce a discrete version of the KiSS 
model, generalizing the stochastic particle model of 
Ref. |24| . As an essential prerequisite, when the num¬ 
ber of particles is large, the model must reproduce 
the FKPP (Fisher-Kolmogorov-Petrovskii-Piskunov) dy¬ 
namics dtO = Dd^6 + r9{l — 9) (i.e., Eq. ([S]) in the dimen¬ 
sional version). The FKPP equation can be derived from 
two coupled reaction-diffusion equations for the concen¬ 
trations 9a and 9b of the species undergoing the auto- 
catalytic reaction 


A + B 


r 



2B. 


( 6 ) 


Summing the dynamical equation for 9a {dt9A = 
DAd^9A — r9A9B) with that for 9b {dt9B = DBd‘^9B + 
r9B9B) gives the FKPP equation, provided Da = Db = 
D and 9a + 9b = 1 (i.e. local mass conservation). We 
discretize these reaction-diffusion equations passing from 
the concentrations 9a.b to a particle description. Thus, 
we consider Na particles of type A and Nb of type B, 
with N = Na + Nb fixed. All the particles undergo 
independent diffusive processes with the same diffusion 
coefficient, D, within the favorable patch [0,L]. Their 
positions Xa (a = l,...,iV) diffuse with diffusion co¬ 
efficient D. As for the interaction between particles, a 
simple way to implement the interaction term r9A9B is 
to impose that a particle of type A changes into B with 
a probability rate Wab depending on the intrinsic rate, 
r, and on the number of B particles within an interac¬ 
tion distance R. This is mathematically formalized by 
the expression 


Wab 


nB{Xa]R) _ nB{Xa]R) 

nav{R) ^ “^Rp 


(7) 


where nB{xa]R) denotes the number of H-particles 
within a distance R from the given A-particle located at 
Xa, and UaviR) = ‘i.Rp the average number of particles 
(regardless their type) within the interaction interval of 
length 2R {p = N/L being the particle density). As dis¬ 
cussed in Ref. [2^ , whenever the local particle density is 
large enough, the probabilistic rule above described con¬ 
verges to the FKPP equation. Finally, we have to fix the 
behavior of A- and B- particles at the boundaries. In 
order to reproduce the boundary conditions of the KiSS 


model and locally ensuring mass conservation, we impose 
the following rules: A particles hitting the boundaries are 
reflected (nutrients are present only within the patch), if 
a B particle hits the boundary it is absorbed (as it can¬ 
not survive outside the favorable patch) and it is replaced 
in-place by an A-particle. 

The A particles are used as “auxiliary” particles in or¬ 
der to recover as a limit the continuous FKPP model. 
From an ecological perspective, considering the A parti¬ 
cles as nutrient, it would be interesting to study the case 
in which A and B particles diffuse with different diffusion 
coefficients. Clearly, in this case the FKPP cannot be ex¬ 
pected to be the proper continuous description. Here we 
limit our discussion to the case of same diffusion coeffi¬ 
cients for the two species of particles. In the Appendinx, 
however, as an example we briefly consider the case of 
non-diffusing A particles. 

From the algorithmic point of view, the above model 
consists of two basic steps: the diffusive step - in which 
all particle positions Xa (a = 1, ..., N) evolve according 
to the dynamics dxa/dt = V 2D rja, with {rja} indepen¬ 
dent normal random variables; and the reaction step - in 
which, for each A-particle, one counts all the H-particles 
within a distance R and changes Am B with probability 
WAB^t- Whenever the time step At is small enough the 
order in which the two steps are performed is irrelevant. 

The main difference between the continuous model and 
the discrete one is that while the (continuous) stationary 
state 9 = 0 becomes unstable above the critical patch 
size, the (discrete) absorbing state Nb = 0 can always be 
reached due to demographic stochasticity. In the follow¬ 
ing sections we will investigate in depth this issue. Here, 
we simply observe that the presence of an absorbing state 
forces us to perform ensemble averages in numerical sim¬ 
ulations. Therefore, our simulations consist of a large 
number, Ur, of statistically independent realizations. As 
for initial conditions, each realization starts with N/2 
particles of each type uniformly distributed within the 
patch. Each realization is then followed until the absorb¬ 
ing state {Nb = 0) is reached. The interaction distance, 
R, is chosen to be small compared to the size of the favor¬ 
able patch but large enough to contain some particles on 
average {Rp > 1), which provides a lower bound on the 
minimum number of particles that one can consider, typ¬ 
ically N > 20. In all simulations we set R = 0.1, D = 1 
and r = 1 so that we have as reference in the continuum 
limit the non-dimensional equation ([S]). Finally, the time 
step. At, is chosen small enough such that the diffusive 
step is much smaller that the interaction distance (i.e. 
V 2D At <C R). Here, we have chosen At = 10 

We notice that the model is robust with respect to 
small changes of R, while if R becomes comparable with 
the patch size relevant differences arise, as discussed in 
the Appendix. There, we also briefly consider the subtle 
effects caused by modifications to the rules specifying the 
dynamics of the individual-based model. 

As an example, we show in Fig.[T]the good agreement, 
at least for large L, between the stationary population 





4 



X 


FIG. 1: (Color online) Density field d(x) at stationarity for 
the continuous KiSS model (black solid curve) compared with 
the local density of B-particles, averaged over = 10 real¬ 
izations, in the discrete model (empty circles). The main 
parameters are: N = 10®, R — 0.1, L = 6 . The (blue) dashed 
curve represents the local density of B-particles obtained with 
a different model with fixed A particles (see Appendix for de¬ 
tails) . 

density, 0(x,t —>■ oo), of the continuous KiSS model and 
the long time (but finite, see Sect. IIIII) particle density, 
obtained by ensemble averaging the spatial distribution 
of B particles in the discrete model. 

III. THE QUASI-STATIONARY STATE 

In the discrete model, the presence of an absorbing 
state (Nb =0) makes population extinction possible for 
any finite N, even when L > Lc- Nevertheless, the time 
taken by the system to be absorbed can be very long 
(for N and L large enough) and, in that case, the system 
reaches a quasi-stationary state. The quasi-stationary 
distribution is defined as the probability distribution con¬ 
ditioned on non-extinction [2^ [H, [H, [sl, HI, [s^ . In gen¬ 
eral, it provides interesting information about the prop¬ 
erties of systems with an absorbing state [2^ [2^ but, 
except in some rare cases, it cannot be evaluated explic¬ 
itly. 

As shown in Fig. [TJ the quasi-stationary distribution, 
apart from fluctuations, essentially reproduces the results 
obtained with the continuous model, at least when L and 
N are sufficiently large. It is then natural to study un¬ 
der which conditions this correspondence holds. In par¬ 
ticular, in this section, we focus on the behavior of the 
biomass per unit length defined as (see Refs, (dol l4ll|) 

Bc{t) = Y [ ( 8 ) 

^ Jo 

The label C simply recalls that this quantity is defined in 
the continuum limit, and allows to distinguish it from the 
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FIG. 2: (Golor online) Survival probability and biomass as a 
function of time for two values of the patch size above the crit¬ 
ical value Lc = tt: L = 4 (a) and L = 3.15 (b), with N = 400 
particles. Top panels show the survival probability Ps{t), see 
text. Bottom panels show the time evolution of the continu¬ 
ous (Beit)) and discrete {B{t)) biomasses, the error bars on 
the latter display the standard deviation, asit). The (blue) 
dotted curves denote the minimum and maximum values of 
NB{t)/N observed in all Ur = 1000 realizations. The inset in 
(b) zooms the framed area in the main figure to better show 
that B{t) > at long times. 

analogous quantity in the discrete case, for which no la¬ 
bels will be used. We observe that Bc{t) asymptotically 
vanishes when L < Lc, due to population extinction, and 
approaches a positive steady value, indicating survival, 
when L > Lc- The asymptotic value = limt_^oo Bc{t) 
can be used as an order parameter to define the extinc¬ 
tion/survival transition. 

In the discrete model, the biomass per unit length (|8]) 
is nothing but the ratio between the total number of B- 
particles and the total number of particles, NB{t)/N. 
Since Ag is a time-fluctuating quantity, we consider its 
ensemble average 



where the brackets (• • ■) denote averaging over different 
realizations with the same initial conditions and condi¬ 
tioning on survival (i.e. at each time t only the realiza- 
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tions with Nsit) > 1 are considered). When condition¬ 
ing on non-extinction, at long times, B{t) approaches a 
definite value (Fig. [2]) which is the average over the quasi¬ 
stationary distribution. 

In the bottom panels of Figure [2^,b we compare the 
time evolution of Be it) and B{t) for two different values 
of L > Lc- For the discrete case we also show, as error 
bars, the standard deviation of the ratio NB{t)/N, 

= ( 10 ) 

In the top panels, we show the probability, Psit) (mea¬ 
sured as the fraction of realizations such that Nb (t) > 1 
at time t), that the population survived up to time t. 

For any given N, when L is large enough (Fig. [2Ji), 
B{t) essentially coincides with the continuum-limit and 
all realizations survive (P^ (t) = 1 in the observation win¬ 
dow) for long time. In these circumstances, the minimal 
values of B{t) (lower blue dotted curve in Fig. [2^) are 
well above zero, implying that very unlikely fluctuations 
are needed to reach the absorbing state. Conversely, as L 
approaches the critical value (Fig. Or) this is no longer 
true. As shown in the top panel of Fig. Or, Ps{t) expo¬ 
nentially decays to zero at a relatively short time. More 
importantly, the agreement between B{t) and the con¬ 
tinuum limiting value gets poorer. In particular, at long 
times, B{t) approaches a quasi-stationary value which is 
larger than the corresponding continuum one (Fig. O^, 
inset). This is clearly a consequence of conditioning on 
non-extinction which leads to an overestimation of the 
biomass, for these values of the patch size. 

The latter observation becomes more quantitative 
when looking at Fig. [3] where we compare, as a function of 
L — Lc, the limiting value B'q with the quasi-stationary 
mean value Bqs, estimated by averaging B{t) over the 
time interval in which it fluctuates around a constant 
value (see Fig. EJ. We have also computed Bqs using the 
algorithm proposed in Ref. [i^ , which directly probes the 
quasi-stationary distribution, obtaining (not shown) in¬ 
distinguishable results. The plateau Bqs ~ const > B^, 
observed when L —^ Lc, is essentially due to conditioning 
on non-extinction. The effect is stronger the smaller is 
the number of particles. The dependence of Bqs on the 
number of particles N clearly demonstrates the depar¬ 
ture from the continuum limit when L approaches the 
critical value. However, for any L > Lc there exists a 
value of N (which diverges as P —^ Lc) above which Bqs 
tends to the continuous value B^. 

When the population survive, it is possibile to show, 
by using a meanfield-like argument, that the biomass of 
the continuous field can be approximated by 

B^iL) « Ai, (II) 

as confirmed in Fig. [31 The idea is to approximate the 
right hand side of Eq. m as d^9 + 9{1 — 0), where 0 de¬ 
notes the steady state spatially averaged value. The first 
eigenvalue of this problem, , is linked to the eigenvalue 
of the linear KiSS model (jl]) via A^ = Ai — 6. To have 



FIG. 3: (Color online) Limiting biomass B'q and average 
biomass in the quasi-stationary state, Bqs, as a function of 
L — Lc for different values of N as in legend. The solid curve 
displays 0.75Ai that corresponds to the behaviour Eq. (EJ. 
The factor 0.75 cannot be catched by the argument exposed 
in the main text. Inset: standard deviation of the biomass 
(O in the quasi-stationary regime, obqs , as a function of N 
for two values (L = 3.5 and L = 7) of the patch size. The 
number of realizations used for the averages ranges from 500 
to 10000 depending on the value of N. 



FIG. 4: (Color online) Quasi-stationary distribution for L — 
Lc = 0.3 and three values of N as in label. The solid curves 
display the Gaussian distribution obtained using the mean 
and standard deviation values measured from data. We ob¬ 
serve a fairly good agreement at least for N large enough. 
When L —>■ Lc and/or N is not large enough the agreement 
ceases to be so good, as one can observe from the left tail at 
N = 500. 


a stationary solution one has to impose A* = 0, so that 
6 — Xi. But at stationarity 0 = B^, which implies dill) . 
Of course this kind of argument cannot be expected to 
work when 9 is close to 1. It is worth observing that, for 
L — Lc <C 1, Eq. (HU implies that the continuous biomass 
is linear in L — Lc- 
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The inset of Fig. [3] shows the standard deviation of the 
biomass in the quasi-stationary state, at varying 

N, for two values of L. At least for L and N large enough, 
the numerical results indicate that ^ We 

notice that the inverse square-root dependence of obqs 
on N implies for the standard deviation of the number 
of i?-particles un = Ngbqs suggesting that the 

distribution of Nb particles is Gaussian, as confirmed in 
Fig. m This observation will be exploited in the next 
section. 


IV. CRITICAL PATCH SIZE AND 
EXTINCTION TIMES IN THE DISCRETE 
MODEL 

In the continuous KiSS model, the critical patch size 
Lc discriminates between populations able to persist 
{L > Lc) and those doomed to extinction (L < Lc). 
For its ecological importance, many efforts have been de¬ 
voted to understanding how this critical length is mod¬ 
ified in less idealized situations, e.g. taking into ac¬ 
count more complex heterogeneous environments 1 ^ 
and/or t he p resence of external factors like fluid advec- 
tion [ 2 ^, IdlL 1^ 1^ . On the other hand, as previously 
mentioned, in the case of the individual-based model ex¬ 
tinction can happen even when L > Lc- It is then nat¬ 
ural to wonder how the survival/extinction transition of 
the continuous KiSS model translates in the behavior of 
the average time to extinction Tg [s^ . which is linked 
to the survival probability Ps{t) as Te = J^Ps{t) dt (d^ . 
Computations of mean extinction times have been mainly 
carried on in non-spatial models, where the dependence 
of the extinction time on the number of individuals 
N at varying the population growth rate has been ob¬ 
tained using diffusive approximations [l^. These results 
have been influential for their rather simple mathemat¬ 
ical formulation and generality. However, it has been 
recently recognized that the diffusive approximation can 
give wrong answers [ 2 ^ [H, [s^ . 

In what follows we examine the dependence of Tg on 
the population size N in the extinction [L < Lc), critical 
(L = Lc) and persistent [L > Lc) regions of the con¬ 
tinuous KiSS model. The different regions, indeed, are 
characterized by unalike functional dependencies of the 
time to extinction Tc on N as shown in Fig. [S] We ob¬ 
serve that quantitatively Tc depends on the initial condi¬ 
tion, however the dependence on N and L, which is here 
explored, is qualitatively independent of the initial con¬ 
dition, provided the initial population is not too small in 
order to avoid spurious fast extinctions. Therefore, we al¬ 
ways considered the initial population Nsit = 0) = N/2. 
An alternative strategy could be to start from the quasi- 
stationary state. However, when N is small and/or 
L < Lc, as seen in the previous section, the quasista¬ 
tionary state is not so well defined from a physical point 
of view. 

Let us first consider the extinction region {L < Lc). In 


this case, for the continuous model, extinction is signaled 
by the density field exponentially going to zero with rate 
given by the first eigenvalue, Ai = l—{Lc/L)'^. In partic¬ 
ular, at sufficiently long times we have that the biomass 
per unit length behaves as Bc{t) ^ with b some 

constant. In the discrete model, it is natural to assume 
that extinction will take place at a time Tc such that 
Bc{Tc) ~ 1/IV, implying, apart from a constant b' , the 
logarithmic dependence on N 

Tc^^AogN + b', (12) 

|Ai| 

which is consistent with the numerics, at least for N large 
enough (Fig. [S]). We observe that the logarithmic behav¬ 
ior m is also found in non-spatial models [ 2 ^. 

Conversely, for L > Lc, Tc displays an exponential de¬ 
pendence on N (Fig. [5|), this is a quite robust feature 



FIG. 5: (Color online) Mean extinction time, Te, vs popula¬ 
tion size, N, for different patch sizes L — Le + SL. Panel (a) 
from bottom to top 5L = —0.3,—0.1,—0.02, 0, 0.02, 0.1, 0.3, 
encompassing the extinction (empty red symbols) and persis¬ 
tence (blue filled symbols) regions, as well as the critical point 
(black half filled symbols) of the continuous model. The con¬ 
tinuous red lines display the prediction m with b' fitted from 
data; the dashed blue lines correspond to the exponential be¬ 
havior m with a fitted. For L = Lc we found a power-law 
behavior, Te ~ , with best fitting exponent 7 = 0.565 (see 

text). Panels (b) and (c) display a zoom of the main figure, 
in appropriate semilog scale, of both the logarithmic (I12II and 
the exponential m regime, respectively. We emphasize that 
the slope in (b) is not fitted but obtained from the first eigen¬ 
value. Depending on the value of N the number of realizations 
ranges from 500 to 5000. 
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observed in non-spatial models [ 1 ^ [H, [s^ where it is 
derived either with the diffusive approximation or more 
rigorous large--/V approximations (see Ref. [25j | for a com¬ 
pact review). Here, the spatial extension complicates the 
use of those well established techniques. However, we ob¬ 
serve that the exponential behavior for the time to extinc¬ 
tion is consistent, at least for L and/or N large enough, 
with the Gaussian behavior of the quasistationary dis¬ 
tribution observed in the previous section (see Fig. n. 
Indeed, we can assume 


Pqs{Nb) oc -exp - 

O'AT \ 


{Nb - {NB)f 

2a‘jf 


(13) 


and it is reasonable to expect that the extinction time 
will be controlled by the small Nb tail of the quasista¬ 
tionary distribution. Roughly, the idea is that if before 
the extinction a well defined quasistationary state sets 
in, then the extinction will take place when Nb ~ 1 [3 
and the time to reach this state will be proportional to 
the inverse of its probability so that Te ^ ^/Pqs{PIb ~ 
1) Ri I/Pqs^Nb —>■ 0). Hence, using (fT^ we obtain 



where we have used that {Nb) = NBqs, with Bqs in¬ 
dependent of N (as is true for large N, see Fig. ID, and 
On = NajSqs ^ ^is from the inset of Fig. |31 

Numerical findings show that the logarithmic and ex¬ 
ponential behaviors are separated by the power-law Te ^ 
N'^ at the critical patch size Lc of the continuous model. 
In particular, our best fit gives 7 = 0.565, which is not 
far from the analytical estimate Tj, 7 ^ found in non- 
spatial stochastic logistic models [^,[^. However, while 
not large, the difference in the value of the exponent is 
nevertheless clearly measurable and might be due to the 
spatial structure of the system under study. 

We conclude briefly discussing the dependence of Tg 
on L (Fig. | 6 ]). When L < Lc, the extinction time Te is 
mainly controlled by the eigenvalue Ai of the continuous 
KiSS model. Indeed, from Eq. using the expression 
of Ai we have 


Te 




logiV, 


(15) 


which is in fairly good agreement with the numerical data 
when N is large (inset of Fig. [5]). 


V. CONCLUSIONS 

We studied an individual-based reaction-diffusion sys¬ 
tem of ecological interest. We focused on a model of a 
discrete population in a favorable patch, surrounded by 
a deadly environment, which recovers the KiSS model in 
the continuum limit. The presence of an absorbing state 



L 


FIG. 6: (Color online) Mean extinction time Te as a function 
of the patch size L for N = 40 (red circles), N — 100 (orange 
triangles) and N = 1000 (blue squares). In the inset Te / log N 
versus L is compared to the asymptotic prediction (I15II (solid 
black curve). The vertical dashed line marks the critical value 
Lc = IT. The number of realizations used for the averages 
ranges from 500 to 5000 depending on the value of N. 


together with demographic stochasticity makes popula¬ 
tion extinction possible even when the continuous model 
would allow for population survival. In such a case, how¬ 
ever, the system attains a quasi-stationary state that we 
have investigated at varying the system and population 
sizes. We have shown that above the KiSS critical patch 
size the biomass in the quasi-stationary state recovers 
the continuum limit value only if the population is large 
enough. On the other hand, when the population has a 
small number of individuals, the link between biomasses, 
dehned in the continuous and particle models, ceases to 
exist. In particular, when the patch size tends to the crit¬ 
ical value the number of individuals required to recover 
the continuum limit diverges. 

Moreover, we have shown that the transition from ex¬ 
tinction to survival translates, in the discrete model, in 
the transition from a logarithmic to an exponential de¬ 
pendence of the average time to extinction, Te, at vary¬ 
ing the population size, N. At the transition, these 
two behaviors are separated by a power-law dependence 
Te ~ N'^ with an exponent definitively different from the 
analytical prediction obtained in the non-spatial logistic 
model. 

We conclude mentioning that it would be interest¬ 
ing in the future to investigate the effects of popula¬ 
tion discreteness in more complex heterogeneous environ¬ 
ments and possibly in the presence of advection 

[^1^.14^. In such cases, besides the survival/extinction 
transition, we expect that the discreteness of the popu¬ 
lation will impact in nontrivial ways on the spatial prop¬ 
agation properties of the population. 














Appendix A: Effects of changes of the microscopic 
rules of the individual-based model 

In this Appendix we consider how two different modi¬ 
fications in the individual-based model affect the system 
dynamics: the first one regards variations of the interac¬ 
tion distance i? - in order to check the robustness of the 
particle model, the second one regards particles’ motility 
and consists in considering fixed A particles - to investi¬ 
gate changes in the continuum limit. 



FIG. 7: (Color online) Average extinction time vs N for 
L = Lc at varying R as from the legend and considering non- 
motile A-particles. The power law behavior Te ~ y^o.565^ 
identifying the extinction/survival transition, is robust for 
variations of the interaction distance around the reference 
value R = 0.1; for much larger values {R = rr/d) the criti¬ 
cal behavior is considerably altered. The same happens if one 
considers a model in which particles of type A do not move 
(filled diamonds). 


The interaction distance i? is a crucial parameter for 
the individual-based model and it is associated to fluctua¬ 
tions of the number of individuals. A significant variation 
of R should thus have important consequences, through 
its impact on the effective growth rate. In Figure 0 we 
present a comparison between numerical results obtained 
in the reference case R = 0.1, used in this study, and 
those obtained with different values of R. For illustrative 


purposes only, we focus on the average extinction time 
as a function of N by fixing L at the critical value Lc- 
As shown in the figure, the power-law scaling Tg ^ N'^ 
(with 7 ~ 0.565), and hence also the value of the critical 
patch size in the discrete model, is robust if the interac¬ 
tion distance varies within the same order of magnitude 
(R = 0.05 and R = 0.2). 

The situation is different when R becomes comparable 
with the patch size: Tg is no longer described by a power 
law (see the case R = 7r/4 in Fig. [7]) meaning that the 
discrete model is no longer at the critical point. Indeed 
at large N the curve Tg = Tg{N) bends towards a log¬ 
arithmic shape as in Fig. [S]when L < Lg. Actually, a 
similar tendency towards smaller values of Tg character¬ 
izes the large N behavior for R = 0.2, although the effect 
is so weak that it is difficult to be seen in the figure. The 
faster extinction is likely due to the fact that increas¬ 
ing R causes the hostile boundary conditions to be felt 
deeper inside the favorable patch. In other words, while 
on average the number of particles within R gets larger, 
the fraction of B particles does not, due to its depletion 
close to the boundary. 

As for the modifications of the motility of particles, 
we study a case in which A particles are fixed, and the 
only other change in the particle model is that, in or¬ 
der to preserve the total number of individuals when a 
B particle is absorbed at the boundary, a new A par¬ 
ticle is introduced at a random position (uniformly in 
the patch). This model has a biological justification in 
the case in which individuals of type A (the nutrient) 
are non-motile (e.g. a plant species) and individuals of 
type B belong to a motile species feeding on it. The 
new model shows a deep alteration in the critical behav¬ 
ior at L = Lg (see the filled diamond symbols in Fig. 0 
and the continuum limit of such a model is not given by 
Eqs. dlllH) (see also Fig. [T]) because particles of type A 
and type B do not undergo the same diffusive process. 
One can expect that in this case the effective growth 
rate is reduced due to less frequent encounters between 
individuals of different types - A particles close to the 
boundaries have less chances to turn into B particles, 
and extinction times result to be smaller with respect of 
the original particle model. This highlights once more 
that the relation between individual-based and contin¬ 
uous reaction-diffusion models is a delicate one, which 
needs to be carefully taken into account when modeling 
the evolution of biological populations. 
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